Splitting of Andreev levels in a Josephson junction by spin-orbit coupling 
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We consider the effect of spin-orbit coupling on the energy levels of a single-channel Josephson 
junction below the superconducting gap. We investigate quantitatively the level splitting arising 
from the combined effect of spin-orbit coupling and the time-reversal symmetry breaking by the 
phase difference between the superconductors. Using the scattering matrix approach we establish a 
simple connection between the quantum mechanical time delay matrix and the effective Hamiltonian 
for the level splitting. As an application we calculate the distribution of level splittings for an 
ensemble of chaotic Josephson junctions. The distribution falls off as a power law for large splittings, 
unlike the exponentially decaying splitting distribution given by the Wigner surmise - which applies 
for normal chaotic quantum dots with spin-orbit coupling in the case that the time-reversal symmetry 
breaking is due to a magnetic field. 
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I. INTRODUCTION 

A Josephson junction is a weak link between two super- 
conductors with an adjustable phase difference <f>. The 
weak link may be a tunnel barrier or a normal metal. 
Fig- [H shows, for example, a Josephson junction consist- 
ing of a small piece of normal metal (a quantum dot), 
connected to the superconductors by a pair of narrow 
constrictions (quantum point contacts). The excitation 
spectrum below the superconducting gap A consists of 
discrete energies, called Andreev levels. In zero magnetic 
field, the energy levels e n are determined by the normal- 
state transmission eigenvalues T n if A <C fr/r^w, where 
Td w is the dwell time of an electron in the normal region 
(before it is converted into a hole by Andreev reflection 
at the superconductor). The relationship is 1 

e n = Ay/l-T n sin 2 (0/2) + 0(A 2 T dw /h). (1) 

Each level is twofold spin-degenerate (Andreev doublet). 

Recently the effect of spin-orbit coupling on Josephson 
junctions became a subject of investigation^^^. This 
is a subtle effect for the following reason: On the one 
hand, in the absence of magnetic fields the normal-state 
transmission eigenvalues T„ are Kramers degenerate be- 
cause of the time-reversal invariance of the normal sys- 
tem. On the other hand, one would expect a breaking 
of the degeneracy of the Andreev doublets because the 
phase difference between the superconducting contacts 
breaks the time-reversal symmetry of the system. Still, 
to leading order in Ard w /^ the one-to-one relationship 
(U) between e„ and T n ensures that the Andreev levels 
remain degenerate for nonzero </>. As was pointed out by 
Chtchelkatchev and Nazaro\», to see a splitting of the 
Andreev doublets as a result of the combined effect of 
spin-rotation symmetry breaking by spin-orbit coupling 
and time-reversal symmetry breaking by the phase differ- 
ence one has to go beyond the leading order in At^ w /H. 




Figure 1: Sketch of a quantum dot Josephson junction: the 
quantum dot (N) is connected to two superconductors (S) by 
point contacts. Spin-orbit coupling splits the energy levels of 
the system when the superconductors have a nonzero phase 
difference. 

This tunable level splitting was exploited in a proposal 
of Andreev qubits for quantum computation 4 . 

In this work we examine the splitting of the An- 
dreev doublets quantitatively by calculating the first or- 
der correction to the energy levels in the small parameter 
Ardw/S. We concentrate our attention on the case when 
the quantum point contacts support one propagating 
mode each. We give a simple relation between the effec- 
tive Hamiltonian for the level splitting of Chtchelkatchev 
and Nazarov 4 and the Wigner-Smith time delay matrix, 

Q = - !S <f, (2) 

where S is the scattering matrix of the normal system. 
As an application, we calculate how the splittings are 
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distributed for an ensemble of systems where the two su- 
perconductors are connected by a chaotic quantum dot, 
assuming that the spin-orbit coupling in the dot is strong 
enough that the dot Hamiltonian can be modeled as 
a member of the symplectic ensemble of Random Ma- 
trix Theory (RMT) 7 - 8 . The present study in the regime 
A <C fi/r^w complements earlier work2^ in the opposite 
regime A > h/r dw . 

The paper is organized as follows. In Sec. [Til we em- 
ploy the scattering matrix approach for calculating the 
first order correction in Ar dw /fi to the Andreev levels, 
and obtain the effective Hamiltonian for the level split- 
ting in terms of the time delay matrix Q. For simplicity, 
we consider the single- channel case in Sec. |TTJ and give 
the multichannel extension in an Appendix. We apply 
our single-channel formula to a calculation of the split- 
ting distribution for an ensemble of chaotic Josephson 
junctions in Sec. IIIIl We conclude in Sec. IIVI with a 
comparison of the splitting distribution of the Andreev 
doublets and the Wigner surmise of RMT. 



II. SPLITTING HAMILTONIAN AND 
WIGNER-SMITH MATRIX 

For energies below the superconducting gap A the 
Josephson junction supports bound states, with excita- 
tion energies given by the roots of the secular equation 1 



Det [1 - a(e) 2 r* A S c (e)r A S h (e)] = 0, 



(3) 



where 
a = cxp 



-iarccos 



e 1 ^ 2 ! . (4] 
e -^/ 2 l > ,[ ' 



and S c (e) and 5h(e) are the scattering matrices of the 
normal system for electrons and holes. They are related 
as 



S h (e) = TS c {~e)T- 



(5) 



where T = ia 2 K is the time-reversal operator for spin- 
1/2 particles. The matrix a 2 is the second Pauli matrix 
acting on the spin degree of freedom and K is the op- 
erator of complex conjugation. Relation (JSJ) reflects the 
fact that in the normal part the dynamics of the holes is 
governed by the Hamiltonianii 



H h = -THJT' 



(6) 



the negative of the time reversed electron Hamiltonian 

We consider the case when the normal part is time- 
reversal invariant, which imposes the self duality condi- 
tion 5 = a 2 S T a 2 on the scattering matrix. (The super- 
script T refers to matrix transposition.) The elements 
of S e (e) change significantly if e is changed on the scale 
of h/r^w, therefore to leading order in Ar^/fi one can 
neglect the energy dependence of S e (e), and take it at 



the Fermi energy, S e (e) ~ 5 e (0). Making use of the 
self-duality of the scattering matrix, and introducing the 
usual block structure 



5 = 



r t' 
t r' 



the secular equation ([3]) can be simplified toi 
1-4^ -tHsin 2 ( - 



Det 



£ 

A 2 



0. 



(7) 



(8) 



From this equation follows the relation (TTJ) between the 
energies and the transmission eigenvalues. 

The correction of order A 2 Td w /^ comes from consid- 
ering the energy dependence of the scattering matrix to 
first order, 5(e) « 5(0) + (dS/de)e. For simplicity, we re- 
strict ourselves here to the case of two single-channel 
point contacts. (The extension to multichannel point 
contacts is given in App. [AJ) For single-channel point 
contacts the self-duality of the scattering matrix implies 

r = pl 2 , r' = p'l 2 , t' = a 2 t T a 2l t = VTu, (9) 

where p, p' are complex numbers, t 2 is the 2x2 unit 
matrix, 1 > T > and U is a 2 x 2 unitary matrix. 
Writing the energy as £q + Se with 



e = A^l-Tsin 2 ^), 



(10) 



and keeping terms up to linear order in the small quanti- 
ties Se — 0(A 2 T(j w //i) and At^/H, one finds the eigen- 
value equation 



Det 



— {o 2 Q T xl a 2 - Qn) sin(0) 



Se 



(11) 



for the energy correction Se. The matrix Q has the block 
structure 



Qu Q12 
Q21 Q22 



(12) 



inherited from the transmission-reflection block structure 
([7|) of the scattering matrix. 

The second term in the determinant (fTTjl shifts both 
eigenvalues by the same amount fe s hift, while the first, 
manifestly traceless term is responsible for the splitting 
±fe S piit of the doublet. We see that the splitting is de- 
termined by the effective Hamiltonian 



sin(0), 



(13) 



with S a traceless Hermitian 2x2 matrix having ma- 
trix elements of order unity. This is the result of 
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Figure 2: A schematic illustration of the splitting of the An- 
dreev doublet as a function of the phase difference <j> for a 
single-channel Josephson junction with spin-orbit coupling. 
The energies are the sum of a degenerate part eo + 5e s hift 
that is even in <j> and a splitting ±fe a piit that is odd in <f>, as 
explained in the text. The maximal splitting is reached at 
(f) = tt/2. 



Chtchelkatchev and Nazarovi Our analysis gives an ex- 
plicit relationiS, between the matrix £ and the time delay 
matrix Q, 



S = -. (osQfiaa - Qn) • (14) 

4Tdw 

This is the key relation that will allow us, in the next 
section, to calculate the level splitting distribution from 
the known properties of the time delay matrix in a chaotic 
system. 

We conclude this section with a symmetry considera- 
tion. The shift fe s hift is even in (f>, just like the zeroth 
order term eo- In contrast, the splitting fe S plit is odd in 
(f>. This is in accord with the symmetry of the Hamil- 
tonian H that gives the full excitation spectrum of the 
Josephson junction. Under time reversal, in our case of 
a time-reversal invariant normal part, it transforms as 
TH((f))T~ l = H(—<p), therefore, for an eigenstate 'F one 
has 

H{4>)T^{-<t>) = E{-4>)T^{-<t>). K > 

An Andreev doublet is therefore of the form 
{s{4>), £■(—</>)}. The decomposition of e{<j>) into even 
and odd parts in <f> amounts to a decomposition of the 
doublet into a degenerate even part and an odd splitting 
part. The resulting <f> dependence of the doublet is 
shown schematically in Fig. O 



III. SPLITTING DISTRIBUTION IN CHAOTIC 
JOSEPHSON JUNCTIONS 

As an application of our general result (|14[) we calcu- 
late how the level splittings are distributed for an en- 
semble of Josephson junctions where the normal part is 
a chaotic quantum dot. We assume that the spin-orbit 
coupling inside the dot is strong enough that the dot 
Hamiltonian can be modeled as a member of the sym- 
plectic ensemble of RMT, i.e. that the spin-orbit time 
r so is much shorter than Td w - 

The splitting distribution can be obtained from the 
known distribution of the scattering matrb*£, and of the 
dimensionless symmetrized Wigner-Smith matrb*i2, 

Qe = -i — S- 1/2 (dS/de)S- 1/2 . (16) 

The distributions of S and Qe are independent^, which 
makes it advantageous to express Q in terms of S and 
Qe- 

Q= T -^S- x ' 2 Q E S l '\ (17) 
n 

In the single-channel case one has 

(18) 

The rates 7„ are distributed according toi£ 

Pill, 7*) oc 1 71 -72I 4 7i72 exp[-4(7! +72)]. (19) 

The distribution of the phases 4> n is£ 

P(0 1 ,0 2 )cx| e *-e*| 4 . (20) 

The matrices of eigenvectors M\ and Mi are members 
of the group Sp(2) of 4 x 4 unitary symplectic matrices, 
and are uniformly distributed with respect to the Haar 
measure of the group^i^. The Haar measure is given as 

d\i oc y\Detg\Iljdxj, (21) 

in terms of the metric tensor <?, defined by 

Tr (dMdM*) = 9ijdxidxj. (22) 

ij 

Here {xi} is a set of independent variables parameterizing 
the Sp(2) matrix M. 

A convenient choice to parameterize Sp(2) is the de- 
composition 

_ / cos(9) MO) W\ (U \ 
M ~\- sin(^) W cos(9) J \0 V J ' 1 6 > 
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Figure 3: Main plot: Distribution of the maximal splitting of 
the Andreev levels (reached at <f> = n/2) in units of A 2 Td w /ft- 
The smooth curve is the prediction of Random Matrix Theory 
calculated from Eq. ((27}, the histogram is the result of a nu- 
merical simulation using the spin kicked rotator. Inset: com- 
parison of the Andreev doublet splitting distribution (solid 
line) and the Wigner surmise (dashed line). For this com- 
parison, the energies are rescaled such that the mean of the 
distributions is unity. 



where W, U and V are SU(2) matrices, and 6 S [0,7r/2]. 
It is seen that the SU(2)®SU(2) factor corresponding to 
the block-diagonal matrix with U and V cancels from 
the spectral decomposition (fT5|) of Qe and S. Using the 
Euler angle parameterization for SU(2), 



U= 



e -i(^a+Vc/)/2 cos(6>,y/2) -eWv-foV 2 sin(6» c/ /2) 
sin^/2) e *(4>u+TPu)/2 cos (^/2) 

(j>u € [0,2tt], ipu £ [0,4tt], 6 v € [0,tt], 

and similarly for the matrices V, W , one finds that the 
Haar measure on Sp(2) corresponding to the chosen pa- 
rameterization is 

d[i(M) <x sin 3 (9) cos 3 (8)d8 JJ sin(0j uIg.iIOjIc,. 

j=U,V,W 

We define the maximal dimcnsionlcss splitting q of the 
Andreev levels (reached at (f> = tt/2) by the formula 



fespiit = gA A ^ dw sin(0). 
The distribution of q is given by 



(26) 



P(q) = J dn(S)dn(Q E )S(q - v/-Det(S)), 
dfj,{Q E ) = d/z(Mi)(i7id72P(7i,72), 
d/j,(S) = d^(M 2 )d(p 1 dip 2 P((pi,f2)- 



(27) 



(q) = 0.181, V<9 2 ) - <<Z> 2 = 0.152 . 
The splitting distribution near zero behaves as 

P(q)~q 2 (q^O). 
For large splittings we find 

P(q)^q- (i (g-^00). 



(28) 



(29) 



(30) 



In order to check our prediction (|27| for the level 
splitting distribution, we have numerically simulated the 
chaotic quantum dot Josephson junction of Fig. [1] using 
the spin kicked rotato r 13 ! 14 . The spin kicked rotator is a 
dynamical model, from which one can extract scattering 
matrices characteristic of chaotic cavities. These scatter- 
ing matrices are given by 



S(e) = V\e' ie - T(\ - V T V)Y Y TV 1 



(31) 



Eq. (|2~T|) can be evaluated numerically. The resulting 
distribution is shown in Fig.[3J The first two moments of 



where T is a 2M x 2M matrix giving the stroboscopic 
time evolution of the model and V is a 4 x 2M projection 
matrix projecting onto the two single-channel point con- 
tacts (the factors of 2 in the dimensions are because of the 
spin). The quasienergy e plays the role of the energy vari- 
able, measured in units of h/to with to the stroboscopic 
time. For a more detailed descri ptio n of this numerical 
model we refer the reader to Ref. YL4 . 

Scattering matrices generated through Eq. (|31[) are in- 
serted into the secular Eq. ([3]), and the roots are found 
by varying the quasienergy. The dwell time in this 
model is Td w = M/2 (again in units of t ). We take 
M = 100 and A = 2 ■ 10~ 4 (in units of h/t ), so that 
Ar dw /fi = 10~ 2 -C 1. By sampling about 10 5 different 
J 7 , P, and 4> we numerically obtain the distribution P(q) 
shown in Fig. [3] together with the analytical result (|27[) . 
The agreement is very good. 



IV. DISCUSSION 

A. Summary 

We have investigated the effect of spin-orbit coupling 
on the subgap spectrum of single-channel Josephson 
junctions. Using the scattering matrix approach and con- 
sidering the energy dependence of the scattering matrix 
to first order we obtained a simple relation, Eq. (|14[) . 
between the effective Hamiltonian governing the level 
splitting and the quantum mechanical time delay ma- 
trix Q = —iS^dS/de. This relation allowed us to find the 
splitting distribution for an ensemble of chaotic Joseph- 
son junctions using the known properties of Q. We ver- 
ified our result numerically by simulating the chaotic 
Josephson junction using the spin kicked rotator, and 
we found excellent agreement. 
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B. Comparison of the splitting distribution with 
the Wigner surmise 

In the inset of Fig. [3] we compare the splitting distri- 
bution of the Andreev doublet with the Wigner surmise 
of RMT4, 



Pw(x) 



32 



x exp 



4a; 2 



(32) 



steps leading to Eq. (fTTj) one arrives at the equation 



Det 



where 



Hn = 



A 2 

2e (0) 



A 2 

H + — K - Se 



-(0)' 



0, 



(Al) 



(A2) 



(For this comparison the energy scale is set such that the 
average splitting is unity.) The motivation behind this 
comparison is the fact that the Wigner surmise is also 
a splitting distribution: as shown in App. [B] it describes 
the distribution of the splittings of Kramers doublets for 
normal chaotic quantum dots with spin-orbit coupling in 
the case that the time-reversal symmetry is broken by a 
magnetic field. 

At small splittings, both P and P\y decay quadrati- 
cally. This quadratic decay is a generic feature of the 
splitting of a Kramers degenerate level due to time- 
reversal symmetry breaking. It follows from the fact that 
the splitting Hamiltonian is a 2 x 2 Hermitian traceless 
matrix without further symmetries and from a power 
counting argument^ similar to the one leading to the 
quadratic decay of Pw 

While at small splittings the two distributions decay 
in the same way, we find qualitative differences in the 
opposite limit. At large splittings P decays like a power 
law in contrast to the exponential decay of P\y [cf. Eqs. 
O and J32])]. 

We attribute the deviation of P from the Wigner sur- 
mise to the nonuniform way in which time-reversal sym- 
metry is broken: While the magnetic field in App. B acts 
uniformly throughout the normal quantum dot, the su- 
perconducting phase difference in the Josephson junction 
acts nonuniformly at the point contacts. 



4°) = AJl-T„sin 2 (^/2), 



(A3) 



and if is a matrix with elements of order T^ w /h. An 
eigenvector of tH with eigenvalue T n is also an eigenvector 
of -Ho with zero eigenvalue. The first order correction to 
the zeroth order energy is the first order perturbative 
correction to this zero eigenvalue. 

We introduce the N x 2 matrices W n and W' n which 
contain the two orthonormal eigenvectors of, respectively, 
t't and t'H', both corresponding to the eigenvalue T n . In 
terms of these matrices we define the matrices q\ n and 
Q2n by 



n = WZQ n W n , q 2n = W^Q 22 W' n 



(A4) 

We find that the shift of the Andreev doublet at el ' is 
given by 



Se 



shift 



A 2 ^ 
4 A 



^1- (ei 0) /A) 2 (Tr q ln + Tr q 2n ) . 



(A5) 

while the splitting fe^ pht is given by the two eigenvalues 
of the traceless Hermitian matrix 

A 2 

H is = — ( CT 2 q\ n o- 2 - gin) sin(0). (A6) 



ACKNOWLEDGMENTS 



Appendix B: SPLITTING DISTRIBUTION FOR 
NORMAL CHAOTIC QUANTUM DOTS 



This work was supported by the Dutch Science Foun- 
dation NWO/FOM. We also acknowledge support by the 
European Community's Marie Curie Research Training 
Network under contract MRTN-CT-2003-504574, Funda- 
mentals of Nanoelectronics. 



Appendix A: SPLITTING HAMILTONIAN FOR 
MULTICHANNEL JOSEPHSON JUNCTIONS 



We generalize the relation (|14[) between the splitting 
Hamiltonian and the time delay matrix to the case that 
each of the two point contacts supports N/2 propagating 
modes. (The single-channel case of Sec. [TT] therefore cor- 
responds to N = 2.) In the multichannel case, after the 



We calculate the splitting distribution of a Kramers 
degenerate level for normal chaotic quantum dots with 
spin-orbit coupling, in the case that the time-reversal 
symmetry is broken by a magnetic field. 

The Hamiltonian of the system is decomposed into two 
parts, 



H = H +A, Hl=Ho, A* =A, 



(Bl) 



where Hq and A are 2M x 2M matrices (the factor of 
two is due to the spin). They satisfy 



THqT- 1 = H , TAT- 1 



-A. 



(B2) 



The matrix Ho models the time-reversal invariant part 
of the Hamiltonian and A is a time-reversal symmetry 
breaking term. 
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The eigenvalues of Ho are doubly degenerate (Kramers 
degeneracy). Considering a doublet with energy E , with 
corresponding eigenvectors u±, u% = Tu\, 

H ui = E Q ui, H a u 2 = E u 2l (B3) 

and treating A as a perturbation, first order degenerate 
perturbation theory leads to the splitting of the Kramers 
doublet by an amount ±(5e sp ii t . We find 

<5e sp iit = ^(u^Au^ 2 + \{ Ul ,Au 2 )\ 2 . (B4) 

For chaotic billiards, the splitting distribution is given 
bj£ 

P(A) = J dU p(U) J dA P(A)S(\ - 5e Bpm ), (B5) 

where U is the matrix of eigenvectors of Hq, distributed 
according to p{U). (The form of p(U) is not needed for 



the derivation.) The matrix A has distribution 

P(A) cx exp (-w 2 Tr A 2 ) , (B6) 

where v is a positive number. Using the fact that P(A)dA 
is invariant under a unitary transformation with the ma- 
trix of eigenvectors of Hq , one finds 



P(X) = J da db dc P(a, b, c)6(\-^a 2 + b 2 + c 2 ), (B7) 
where 

P(a, b, c) cx exp[-2v 2 (a 2 + b 2 + c 2 )]. (B8) 

After changing to polar coordinates the integral (|B7[) can 
be evaluated straightforwardly, and after rescaling from 
A to x, defined by J dx P(x)x = 1, one arrives at the 
Wigner surmise 
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